function [momentsScaled,nameMomentsScaled] = getScaledMoments(data,autoLagsIdx,inclMoms)

% Number of observations (T) and number of observables (ny)
[T,ny] = size(data);

% The moments
maxLag   = max(autoLagsIdx);
autoEyy  = zeros(ny,ny,maxLag);
numLags  = size(autoLagsIdx,2);
Ey       = sum(data,1)'/T;
Eyy      = corr(data);
stdy     = std(data,0,1);
for i=1:ny
    Eyy(i,i) = stdy(i);
end
for i=1:maxLag
    autoEyy(:,:,i) = corr(data(1:T-i,:),data(1+i:T,:));
end

% We collect the moments
momentsScaled = collectMoments(Ey,Eyy,autoEyy,autoLagsIdx);

% Removing potential moments
momentsScaled     = momentsScaled(inclMoms==1);


if nargout > 1
    % We construct a cell with general labels of the moments
    for i=1:ny
        nameMomentsScaled{i,1} = ['E[y_t(',num2str(i),')]'];
    end
    idxMom = ny;
    for i=1:ny
        for j=i:ny
            idxMom = idxMom + 1;
            if i == j
                nameMomentsScaled{idxMom,1} = ['std[y_t(',num2str(i),')]'];
            else
                nameMomentsScaled{idxMom,1} = ['corr[y_t(',num2str(i),')y_t(',num2str(j),')]'];
            end
        end
    end
    for k=1:numLags
        for i=1:ny
            idxMom = idxMom + 1;
            nameMomentsScaled{idxMom,1} = ['corr[y_t-',num2str(autoLagsIdx(1,k)),'(',num2str(i),')^2]'];
        end
    end
    % Removing potential moments
    nameMomentsScaled = nameMomentsScaled(inclMoms==1);
end





end